function [ Alpha, Beta ] = alphaBeta( S, r, sigma )
%ALPHABETA Summary of this function goes here
%   Detailed explanation goes here

% Now get the alpha and beta vector
S_i = S(2:end-1);
S_if = S(3:end);
S_ib = S(1:end-2);
% Central difference
A_c = (sigma^2*S_i.^2)./((S_i-S_ib).*(S_if - S_ib))-(r*S_i)./(S_if-S_ib);
B_c = (sigma^2*S_i.^2)./((S_if-S_i).*(S_if - S_ib))+(r*S_i)./(S_if-S_ib);

% Forward difference
A_f = (sigma^2*S_i.^2)./((S_i-S_ib).*(S_if - S_ib));
B_f = (sigma^2*S_i.^2)./((S_i-S_ib).*(S_if - S_ib))+(r*S_i)./(S_if-S_i);

% Backward difference
A_b = (sigma^2*S_i.^2)./((S_i-S_ib).*(S_if - S_ib))-(r*S_i)./(S_i-S_ib);
B_b = (sigma^2*S_i.^2)./((S_if-S_i).*(S_if - S_ib));

% Determine central or upstream discretization to use
if (sum(A_c<0) == 0 && sum(B_c<0) == 0)
    Alpha = A_c;
    Beta = B_c;
else
    % Find where A_c is negative and change both node to upstream
    negA_c = A_c < 0;
    % Get the upstream values
    if(sum(A_f<0) == 0 && sum(B_f<0) == 0)
        A_u = A_f;
        B_u = B_f;
    else
        negA_f = A_f < 0;
        A_u = A_b.*negA_f + A_f.*(1-negA_f);
        B_u = B_b.*negA_f + B_f.*(1-negA_f);
    end
    % Now correct Alpha and Beta
    Alpha = A_u.*negA_c + A_c.*(1 - negA_c);
    Beta = B_u.*negA_c + B_c.*(1 - negA_c);

end
    
end

